Computational exploration of novel ROCK2 inhibitors for cardiovascular disease management; insights from high-throughput virtual screening, molecular docking, DFT and MD simulation

Cardiovascular disorders are the world’s major cause of death nowadays. To treat cardiovascular diseases especially coronary artery diseases and hypertension, researchers found potential ROCK2 (Rho-associated coiled-coil-containing protein kinase 2) target due to its substantial role in NO-cGMP and RhoA/ROCK pathway. Available drugs for ROCK2 are less effective and some of them depict side effects. Therefore, a set of novel compounds were screened that can potentially inhibit the activity of ROCK2 and help to treat cardiovascular diseases by employing In-silico techniques. In this study, we undertook ligand based virtual screening of 50 million compound’s library, to that purpose shape and features (contain functional groups) based pharmacophore query was modelled and validated by Area Under Curve graph (AUC). 2000 best hits were screened for Lipinski’s rule of 5 compliance. Subsequently, these selected compounds were docked into the binding site of ROCK2 to gain insights into the interactions between hit compounds and the target protein. Based on binding affinity and RMSD scores, a final cohort of 15 compounds were chosen which were further refined by pharmacokinetics, ADMET and bioactivity scores. 2 potential hits were screened using density functional theory, revealing remarkable biological and chemical activity. Potential inhibitors (F847-0007 and 9543495) underwent rigorous examination through MD Simulations and MMGBSA analysis, elucidating their stability and strong binding affinities. Results of current study unveil the potential of identified novel hits as promising lead compounds for ROCK2 associated with cardiovascular diseases. These findings will further investigate via In-vitro and In-vivo studies to develop novel druglike molecules against ROCK2.


Introduction
Cardiovascular disorders (CVDs) include heart and blood vessels related diseases such as cerebrovascular diseases, coronary heart diseases, hypertension, stroke, arrhythmia and other conditions caused by high blood pressure and cholesterol, diabetes mellitus, physical inactivity, alcohol consumption, smoking, and obesity etc. Aging, stress, microRNAs, and diet also have its own role in CVDs as it causes several changes in human heart [1,2].Risk of developing CVDs is three times higher if one's parents had heart diseases and it is most commonly found in men than in pre-menopausal women [3].Cardiovascular disorders are the leading cause of mortality worldwide and most deaths are due to heart attacks and strokes.In western countries the situation is worse, about sixty-two million people died with the cardiovascular diseases and fifty million people with hypertension.The situation is going to be worse in Pakistan also.According to WHO report 2019, worldwide mortality rate due to CVDs is 32% [4].
Cardiovascular diseases are associated with the dysfunction of NO-cGMP (Nitric oxidecyclic 3'-5' guanosine monophosphate) pathway (Fig 1) and this pathway act as target for treating different cardiovascular diseases.NO (Nitric oxide) act as relaxing factor due to which recognized as vasodilator, causing smooth muscle relaxation [5].Calcium calmodulin complex leads to activation of eNOS (endothelial nitric oxide synthase) which oxidizes the L-arginine into L-citrulline to produce NO as well as it can be produced nonenzymatically.NO targets soluble guanylyl cyclase (sGC) by binding to its heam moiety which is in smooth muscle cells and aids in the relaxation of surrounding smooth muscle cells.Activated sGC increases the concentration of cGMP in tissues which induces the relaxation of vascular tone and in this way there is increase in blood flow which lowers the blood pressure [6].cGMP the secondary messenger, play major role in cardiovascular system even minor disturbance in any step leads to malfunctioning of cardiovascular system such as hypertension, cardiac hypertrophy, myocardial infarction as well as the heart failure [7].cGMP also known as cGK (cGMP dependent protein kinase) that found in wide range of eukaryotes [8].cGKI is most often expressed cGK in cardiovascular systems while cGKII not expressed in cardiac and vascular myocytes [9].There is calcium dependent and independent contraction and relaxation in cGMP/cGKI pathway when activated MLCK phosphorylate MLC and MLCP dephosphorylate MLC respectively [6].
For the treatment of cardiovascular diseases such as hypertension, cerebral ischemia, vascular inflammation, arteriosclerosis, atherosclerosis, coronary vasospasm, and stroke, ROCK2 was used as therapeutic drug target because ROCK2 inhibit myosin binding subunit (MBS) of MLCP to regulate MLC's phosphorylation at Ser19 residue instead of dephosphorylation and cause calcium independent contraction instead of relaxation.Therefore, both MLCK and MLCP phenomena promote contraction and lead towards different CVDs (Fig 2).If ROCK is inhibited, MLCP promote dephosphorylation of MLC and cause relaxation and in this way rate of hypertension will be reduced [10].ROCK2 relate to coronary artery diseases and hypertension as well as it plays vital function in the contraction of smooth muscle cells than ROCK1 [11].Genes of both ROCK1 and ROCK2 have different subcellular localizations where they expressed.Except brain and muscles, ROCK1 is widely expressed in the liver, kidney, lung, spleen, and testis.While ROCK2 is abundantly expressed mostly in brain, heart, muscles, and placenta.All these results confirm the fundamental role of ROCK2 in cardiovascular diseases [12].So, a significant target ROCK2 was employed for the treatment of numerous cardiovascular diseases i.e., hypercontraction, atherosclerosis, hypertension, and heart failure [13].
ROCK2 (Rho-associated coiled-coil-containing protein kinase) is a serine/threonine protein kinase that is downstream effector of RhoA.ROCK2 also known as p164 ROCK-2 or ROKa, ROCK-II and consist upon 1388 amino acids [14].It is positioned on chromosome 2 and consists upon three major domains; at the N-terminus of ROCK2 there is a catalytic kinase  [16].This conformational alteration leads to kinase activation [13].
Fasudil, Y27632 and ripasudil are significant ROCK inhibitors that have been endorsed for the treatment of hypertension, cerebral vasospasm, and glaucoma respectively.In addition, there are Indazole and amidoindazole available as ROCK inhibitors to treat breast cancer.DW1865 is another ROCK inhibitor which has shown good cell potency in fiber formation.Dual ROCK1 and ROCK2 inhibitors have been linked to problematic side effects as well as some investigators also hypothesized that ROCK inhibitors would be less effective for the treatment of cardiovascular diseases [17].Because of these issues it is necessary to predict novel compounds with stronger binding affinity for ROCK2 to treat CVDs.High selectivity for ROCK2 could result into a product candidate with improved tolerability which could allow for long-term systemic use [18].

Material and methods
The graphical workflow of adopted methodology was illustrated in Fig 4.

Structure refinement and evaluation
X-ray crystallographic structure of targeted receptor was fetched from RCSB PDB depicting 0.281 R-value free and 2.79 Å resolution while amino acid sequence of target protein was retrieved from UniProt database in FASTA format and BlastP was run to determine phylogeny of query protein (ROCK2).The Expected value, identity, similarity, and query coverage values help to identify evolutionary relationship between different sequences as well as members of gene family.Expasy ProtParam tool [19] was used to analyze target protein's physiochemical properties such as molecular weight, theoretical pI, atomic composition, instability index, GRAVY (grand average of hydropathicity), and amino acid composition.Domains are accountable for a specific function, which contributes to the overall role of a protein.InterPro database analyses protein sequences functionally by categorizing them into families and anticipates the existance of domains, repeats, binding, and active sites [20].Prediction of protein subcellular localization is an essential step to understand protein function.To anticipate subcellular localization of eukaryotic proteins, LocTree3, DeepLoc-1.0,ESLpred, and HSLPred tools were used which worked on deep neural networks relying solely on their amino acid composition, dipeptide composition, and physicochemical properties.CASTp [21] was used for predicting active site residues that were further validated by site finder function.
Missing residues were added by utilizing UCSF Chimera's built-in discrete optimized protein energy (DOPE) loop modeling protocol and evaluated modeled structure on the bases of zDOPE and GA341 score.GalaxyRefine, and 3Drefine were employed for refinement of loop and terminus regions using ab initio modeling combined with atomic-level energy minimization.SAVES Server (ERRAT, Verify3D, PROVE, and PROCHECK) was used to validate and examine the stereochemical quality of protein structure.ERRAT and VERIFY 3D distinguishes between correctly and incorrectly folded regions of protein structure [22,23].PROCHECK performs a thorough examination of protein structure's stereochemistry as well as draws attention to regions of protein that may require further investigation while PROVE checks the deviation using standard parameters [24].

Pharmacophore query generation and validation
55 co-crystal structures of ROCK2 were fetched from PDB (S3 Table in S1 File), and 18 active compounds were selected based on IC 50 values and their physiochemical properties.Common features were found to build pharmacophore query that depict bioactivity towards target protein.Between two pharmacophoric features, a maximum 10 Å distance was allowed.Validation of the pharmacophore model was accomplished via Receiver Operating Characteristic curve (ROC) yielded area under curve (AUC) and enrichment factor.For this purpose, 18 cocrystal compounds were used as an active control, and a decoy dataset of 619 compounds from ChEMBL database was utilized as a negative control.Physiochemical properties were set as a threshold to filter chemically distinct compounds from ChEMBL and generated 10 conformers of each decoy molecule using Avogadro [25].Target protein docked with active and decoy molecules and docking scores were further used to generate ROC curve.AUC access test's overall accuracy and values ranges from 0-1, where 0 shows perfectly inaccurate test and 1 denotes thoroughly accurate test.An AUC of 0.5 (black diagonal line) indicates no differentiation, 0.7-0.8deemed acceptable, 0.8-0.9considered excellent, and values beyond 0.9 as remarkable.

Pharmacophore-based virtual screening of prepared library
50 million drug-like and lead-like chemical compounds were retrieved from different databases like Asinex, Binding Database, Chembridge, ChemDiv, Drug Bank (FDA-approved molecules), Maybridge, and Zinc Database.Pharmacophore query employed to screen prepared libraries via virtual screening protocol.Virtual screening was carried out to locate hit compounds that share chemical and structural similarities with pharmacophore model and could inhibit target protein's activity.

Molecular docking
Screened compounds were subjected to molecular docking in the active site of ROCK2 to find out best binding pose under physiological conditions (300K, PH 7, 1 atm, 0.1 salt concentration), and for this purpose detailed conformational analysis was undertaken employing multiple docking algorithm and molecular graphics program were used for good quality 2D visualizations.It provides results on the basis of empirical and force field-based scoring functions.Both receptor and hit compounds were prepared by protonation using Protonate3D module and energy minimization by using Amber ff10 and EHT (Extended Hu ¨ckel Theory) force field for macromolecules and ligands respectively.AMI-BCC charges are projected for small molecules as well adjust hydrogens and lone pairs as required.Both protein and hit compounds were saved into their respective file formats then site-specific docking was processed by implementing triangle matcher protocol with rigid receptor docking method.London dG the initial scoring method evaluates initial 30 ligand placements while GBVI/WSA dG (generalized-born volume integral/weighted surface area) scoring function refines best 5 poses.The docked poses of ligands were further evaluated via binding mode, S score, and RMSD values.PyMOL was used for 3D visualization of receptor-ligand while Discovery Studio for 2D depiction.

Evaluation of pharmacokinetic properties
A successful lead needs to be more druglike for successful drug discovery so top compounds were screened based on drug-likeness, ADME (Absorption, Distribution, Metabolism, and Excretion), and toxicity characteristics through PkCSM, ADMETlab, and SwissADME which calculate important drug like descriptors as well as used for predicting mutagenicity, immunotoxicity, hepatotoxicity, and carcinogenicity.The druglike compounds should have to follow certain rules such as Lipinski's, Veber's, Ghose's, Egan's, and Muegge's rules.These rules trace druglike molecules on the basis of 2D depiction (physical, chemical properties) and molecules must have to fulfill up to 2 of these rules otherwise represent some ADMET obstacles.Drug likeness is an important concept in drug development as it helps in identifying compounds of good pharmacokinetic properties [26].These approaches help to reduce the number of compounds that need to be tested in clinical trials as well as save time and resources in drug development.After ADMET and drug-likeness analysis, we will have compounds whose chances of selection in clinical trials increases.

Bioavailability and bioactivity score prediction
Bioactivity and bioavailability scores are significant factors in drug discovery/development and depict overall potential of lead compounds.These factors help researchers to prioritize and optimize compounds for further stages of evaluation.Bioavailability determines compound's effectiveness via physicochemical properties, such as solubility and permeability via SwissADME.To predict biological activity of selected compounds, Molinspiration (www.molinspiration.com)software was employed.Drug molecule's bioactivity can be impacted by several parameters, including its chemical structure, potency, and selectivity for a particular target.The drug is expected to bind primarily with target proteins i.e., enzymes, ion channels, and receptors.A drug molecule with high potency and selectivity for a particular target is more likely to deliver a therapeutic benefit with fewer side effects than a molecule with lesser potency or selectivity [27].Bioactivity score of top 6 complexes were predicted against human receptors such as GPCR (G protein-coupled receptors), kinases, ion channels, proteases, and enzyme inhibitors.

DFT analysis
To analyze the molecular geometry, conformation, energetics and reactivity of top selected compounds, density functional theory (DFT) analysis was performed via Gaussian 09 and GaussView 06 package.HOMO and LUMO energies were estimated to calculate the energy gap between the highest occupied and lowest unoccupied molecular orbitals.To identify electrophilic and nucleophilic functional sites of A1 and A2 compounds, MEP (molecular electrostatic potential) was carried out.For this, top compound's geometries were optimized by applying B3LYP method with 6-3IG** basis set.To figure out stability and reactivity of A1 and A2 compounds, HOMO, LUMO, energy gap (E g = ELUMO-E HOMO ), electronegativity, electrophilicity index, electrochemical potential, hardness and softness parameters were computed.

MD simulation
200 ns MD simulation was employed to study protein-ligand interactions at atomic level and conformational changes during binding process by applying Newton's laws of motion via Desmond from Schrodinger suite [28].Both systems (A1 and A2 complexes) were preprocessed via protein preparation wizard module in Maestro using default settings while whole system was set up by "system builder tool".Minimize energy of both systems then equilibrate in the orthorhombic box (10Å×10Å×10Å) of TIP3P water model, neutralized by adding 0.15M NaCl concentration and heating gradually from 200K to 250K and then 300K with 1 atm pressure [29].Receptor-ligand interactions were examined via simulation interaction diagram tool.The MD simulation was successfully conducted using Dell T7810 Precision Workstation with an Intel Xeon E5-2670 v3 processor, 64 GB RAM, and Zotac GeForce RTX 3070 GPU on a system running Ubuntu 20.04.1 LTS.

MMGBSA calculations
MM-GBSA (molecular mechanics generalized Born surface area) analysis was conducted via prime module of Schro ¨dinger to evaluate binding free energy (ΔG bind ) of both systems (ROCK2 complexed with A1 and A2) [30] throughout the simulation time.Counter ions were stripped, VSGB solvent model with OPLS 2005 force field were employed along rotamer search techniques to calculate ΔG bind .The calculations were based on MD trajectory frames taken at 10ns intervals after MD run.Eq (1) The total binding free energy is the difference between the energy of protein-ligand complex (G complex ) and free energy of individual protein (G protein ) and inhibitor (G ligand ).ΔG bind calculated via Eq (2), a combination of molecular mechanics gas phase energy (ΔE gas ), solvation free energy (ΔGs ol ), and entropy (TΔS) [31].
Gas phase interaction energy (ΔE gas ) is equal to van der Waals, and electrostatic energy while internal energy is neglected (Eq (3).(Eq (4) Solvation free energy (ΔGs ol ) is the sum of polar (generalized born) and nonpolar energy (solvent accessible surface area) [32].TΔS denotes conformational entropy equivalent to translational, rotational and vibrational entropy changes upon binding (Eq (5)).

Results and discussion
The accession number, PDB, and UniProt ID of target protein are NP_004841.2,6ED6, and O75116 respectively.ROCK2 has Cyclic-C2 global symmetry while global stoichiometry is Homo 2-mer-A2, and physiochemical properties are listed in S1

Pharmacophore quality assessment via ROC curve
Ligand-based pharmacophore model contains four pharmacophoric features i.e., F1 and F4 for aromatic functional group, one nitrogen as HBD (F2), one oxygen as HBA (F3) on the basis of PCH-All scheme with 1.A validated pharmacophore model was used to perform virtual screening of prepared library of *50M compounds that yielded 5540 compounds.Then these compounds were filtered on the basis of LIPINISKI'S "rule of five" (MW<500, HBD< 5, HBA< 10, logP < 5), and 2085 hits were retained.Further compounds were screened based on hit rate (80%) and a cutoff of RMSD values of more than 1.5Å.Finally, 2000 best hits retrieved will likely bind with the same active site of ROCK2.

Binding poses and interactions
Molecular docking was performed to further screen hit compounds and establish binding modes of selected molecules that may help to treat the disease by inhibiting ROCK2.Before this, docking protocol validated via ROC curve which depicts false positive and true positive fractions on x-axis and y-axis respectively (S7 Fig in S1 File).Based on receiver operating characteristic curve (ROC curve), 0.853 AUC (area under curve) value was observed, and enrichment factor in top 1% was observed 8.88.The results suggest that docking protocol did not produce false-positive results.
Ligand efficiency is a useful parameter to select lead compound and in the optimization process.Table 1 depicts ligand efficiency, RMSD, and different types of interactions ligand made with receptor along their distances.Fig 9 shows comparison of binding affinity and RMSD of first 6 hit compounds along their correlation analysis which depict interdependence of ligand efficiency and drug score.Compounds were further analyzed by superimposing them on reference compound and S6 Fig in S1 File.illustrates superimposed view of reference compound with two top hits and their surface mapping along 2D depiction.On the basis of results, it could be anticipated that both A1 and A2 compounds favorably bind to the active site of ROCK2 and particularly Lys190, Leu195, and Asp206 residues play significant role to inhibit the activity of ROCK2 receptor.

ADME analysis, drug likeness, and toxicity assessment
Poor ADMET profile is responsible for the dropout of lead compounds during clinical trials.So, it's necessary to trace these issues at the early stages of drug discovery.For this, top 6 compounds carried out for ADME and toxicity analysis.All compounds lie within permissible range specified in RO5 (HA, MW, HBD, HBA, and Log P) and no compound violates Veber, Egan, and Muegge's guidelines.As the successful lead must satisfy the ADMET properties.Here A1 and A2 compounds satisfy ADMET properties along synthetic accessibility score of 3.59, 3.47 with high GI absorption and no blood-brain barrier (Table 2).According to toxicity profile analysis, compounds didn't show tumorigenic, mutagenic, and reproductive effects.As well they didn't show skin sensitization and irritation.Toxicity analysis of selected compounds were enlisted in

Bioactivity score prediction
The lead compound with high bioavailability score will be able to reach its target site in adequate amounts to exert therapeutic effect.Bioactivity scores provide details of drug binding for various receptors such as nuclear receptor ligand, protease, and enzyme inhibitors.Table 3 illustrates bioactivity scores, and all compounds lie within -5.0 to 0.0, showing moderately active compounds while some depict scores above zero illustrating the good activity of compounds.The values clearly indicate that lead compounds possess properties required to act as potential drugs.

DFT calculations
Quantum computational methods, frontiers molecular orbitals (HOMO and LUMO) predict stability and chemical reactivity on the basis of highest and lowest excited states (Fig 12A -12D).To understand inhibitory potential of both compounds, energy gap and other descriptors were computed (Table 4).A1 compound has 0.14735 kcal/mol energy gap with 2.005eV hardness and 0.496eV softness while A2 illustrates 0.11685 kcal/mol Eg, 1.59eV hardness and 0.628eV softness.Both compounds have shorter energy gaps and hardness along higher softness values indicating more reactivity and stability.Both compounds depict increased HOMO (donate electron) and decreased LUMO (accept electron) energies reflect the development of stable interaction and significant receptor binding.Both compounds have antioxidant ability as they show 1.87eV and 2.15eV ionization potential for A1 and A2 respectively.Electrochemical potential values for A1 (-3.875eV) and A2 (-3.74eV) revealing their charge transfer capacity to receptor while electrophilicity index classified both compounds as strong electrophiles.Detailed electrophilic and nucleophilic sites illustrated graphically by molecular electrostatic potential (MEP) which shows detailed geometries and favorable points on compounds that influence the establishment of interactions with receptor (Fig 12E and 12F).MEP highlights favorable spots via different color schemes as the negative areas in red color, positive area in blue and zero potential areas with green color.

MD simulation
To visualize protein-ligand stability and structural constancy of both selected (Chem-Div_F847-0007 (A1) and DrugBank_15985904 (A2)) systems and apo protein, MD simulation of 200 ns was carried out.RMSD, RMSF values of both complexes and apo protein were These findings demonstrate that compared to apo-protein, there are no major structural differences and increase in P-RMSD value while A1 system is more stable than A2.
The observation was further supported by RMSF analysis which determines the significant rigid and flexible regions to determine its functionality and compactness.According to Fig 14A1 and 14B1, atom wise RMSD of A1 lies within 1.2-2.8Å while L-RMSF for A2 display higher peaks within the range of 2-5 Å while both hit compounds contain 8 RBs.(Fig 14A2 and 14B2) Residue wise RMSF of A1 maintain 0.8-3 ± 0.5 Å except one higher peak and P-RMSF for A2 was less than 1-5 Å.The higher RMSF areas indicate presence of loop/coil regions while alpha helices/beta-strands make it rigid with lower fluctuations.The green Table 2. ADMET analysis of best selected compounds along some toxicophoric rule's alert.The prediction values transformed into six symbols: 0-0.1(-), 0.1-0.3(-),0.3-0.5(-),0.5-0.7(+),0.7-0.9(++),and 0.9-1.0(+++).5 like RMSD was 3.0 Å and 2.0 Å for both A1 and A2 respectively.To measure the extendedness of both ligands, rGyr for both ligands were calculated which lies within 4-5.0 Å and 5.4-6.3Å respectively.SASA was within the range of 50-200 Å 2 for best hit compound which stabilize at 45ns throughout the simulation but in case of second-best compound it lies within 160-320 Å 2 .MolSA and PSA for A1 were 380-440 and 60-90 Å 2 and in case of A2 it lies between 412-424 and 184-200 Å 2 respectively.

Comparative conformational analysis across selected trajectories
To get insights regarding time-evolution/conformational changes as well as any significant alterations in key structural parameters, conserved and reformed ligand-amino acid bindings,  while at 200ns frame, Ala205 make significant H-bond with hinge region and cause a loss of initial H-bond with Asp206.In system A2, it is notable that residue Arg105 form H-bond occurs exclusively in 200 ns frame, while Asp107 form H-bond across all three frames.Altogether, A1 system fully occupied the active site while A2 compound demonstrates a partial displacement at 100ns frame.

Binding free energy landscape
Binding free energies were determined for MD trajectory frames of both systems which aid in the design and development of medication.According to Fig 17, both systems possess higher negative binding free energy values i.e., (ΔG bind ) = -54.420kcal/mol and -62.272 kcal/mol respectively and follow the cutoff value of -50 kcal/mol which indicate better ligand binding in receptor's active site.According to detailed energy profile, van der walls (-49.606,-49.376 kcal/ mol) contribute more than other energy types while lipo (-17.64,-23.04 kcal/mol), H-bond (-1.95, -1.66 kcal/mol) and coulombic (-9.75, 2.30 kcal/mol) energy types make significant contributions.These results show consistency with earlier findings but before using them as drugs, preclinical and clinical evaluation needed towards ROCK2 receptor.

Conclusion
The current study makes use of the ROCK2 target protein and computational and medicinal chemistry methods to anticipate prospective lead compounds to inhibit ROCK2 and treat cardiovascular diseases.Pharmacophore validation followed by the high throughput virtual screening of *50M small molecules from different databases, top hits docked and revealed similar binding poses which were reported for reference compound.The compounds displayed excellent ADMET profile, HOMO-LUMO energies, bioavailability, and bioactivity  scores.MD simulations of top systems A1-ROCK2 and A2-ROCK2 revealed stability and didn't undergo large conformational changes.MMGBSA studies revealed strong binding affinity and retaining the activity.In-silico approaches help to retrieve potential compounds A1 and A2 within less time and can be used as promising drug molecules after further optimization, preclinical and clinical trials.

Fig 1 .
Fig 1. Human NO-cGMP (Nitric oxide-cyclic 3'-5' guanosine monophosphate) pathway where Ca 2+ binds with calmodulin and makes a complex that activates eNOS to produce NO by oxidizing L-arginine to L-citrulline.NO targets sGC to increase the production of cGMP which starts a cascade.https://doi.org/10.1371/journal.pone.0294511.g001

Fig 3 .
Fig 3. ROCK2 activation mechanism via RhoA independent and dependent ways.In RhoA-independent way, arachidonic acid binds to PH domain or Caspase 2/Granzyme B cleaves in between the RBD and PH domain while in RhoA-dependent way, activated RhoA (RhoA-GTP) binds with the Rho binding domain to activate ROCK2 that further phosphorylate MYPT1 subunit of MLCP and resultantly MLCP become inactive.So, MLCP is unable to dephosphorylate MLC and increase the actomyosin assembly and contractility.https://doi.org/10.1371/journal.pone.0294511.g003

Fig 4 .
Fig 4. Graphical workflow of research.https://doi.org/10.1371/journal.pone.0294511.g004 Fig 5 depicts 3D structure of receptor along the Ramachandran plot and highlights amino acids that are in disallowed regions (white color).According to Ramachandran plot, 95.4% of residues lied in the core regions (red color), 4.3% of residues in additional allowed region (yellow) of torsion values, and 0.3% of residues in generously allowed regions (light yellow) with no amino acid belongs to active site of ROCK2.(S1A Fig in S1 File).Overall quality factor of protein was 93.0748%, 0.26 Å RMSD, and 1.5 MolProbity score indicating physical correctness of all atoms within protein structure.S1B Fig in S1 File.represents correctly and incorrectly folded regions of protein and score ranges from -1 to +1.Approximately, 89.26% amino acids have averaged score > = 0.2 and 0.471 Z-score.S2 Fig in S1 File depicts domains and functional sites of receptor and the main catalytic kinase domain of ROCK2 illustrated in pink color from residue 92 to 354.DeepLoc-1.0,ESLpred, and HSLPred web servers predict that target protein is present in cytoplasm of cell with 0.653 likelihood and 97% accuracy (S3 Fig in S1 File.).Binding pocket residues of receptor identified through CASTp and site finder function are labeled and displayed in Fig 5A.S2Table in S1 File represents binding pocket residues of ROCK2.

Fig 6 .
Fig 5. (A) 3D structure of ROCK2 with labeled active site residues along the (B) Ramachandran Plot exhibiting different regions specified with different colors.Most favored regions with red color, additional allowed regions with yellow color, generously allowed regions with light yellow color, and disallowed regions with white color.https://doi.org/10.1371/journal.pone.0294511.g005

Fig 7
represents top 6 hit compounds while Fig 8 revealed interaction mechanism of topranked poses of 6 hits which bind in the same binding site as the reference compound 3SG.A1 demonstrated best binding score (-14.2884kcal/mol) with 1.09Å RMSD among other 14 compounds.
Fig 10 depicts heatmap of ADME properties of top six compounds ranging from low (blue) to high (yellow).
Fig 11A while Fig 11B depicts distribution and correlation between oral rat acute toxicity and carcinogenicity.According to Fig 11, top A1 and A2 compounds have acceptable toxicity profile, and can be used as lead compounds for further analysis.

Fig 9 .
Fig 9. (A) Comparison of Binding affinity (ΔG) and RMSD of best six hits with the ROCK2 (B) Correlation analysis between Ligand efficiency (LE) and Drug score of best hits via Python Seaborn module.https://doi.org/10.1371/journal.pone.0294511.g009 vertical bars display strong intermolecular ligand-receptor interactions with lowest fluctuations which depict high stability at the active site.Fig 14C1 illustrate RMSF for apo-protein which lies within *4 Å while Fig 14C2 depicts 3D structure of target protein along 3 highlighted flexible regions (R1-R3).It can be observed that there is no major shift recorded in apo and selected system's RMSF, and ROCK2 maintain compactness over 200 ns simulation.Secondary structure elements monitored throughout the simulation.(S8 Fig in S1 File).

Fig 12 .
Fig 12. Electron density regions HOMO and LUMO (A, B) A1 compound and (C, D) A2 compound.E and F display molecular electrostatic potential of A1 and A2 compounds respectively where dark red color represents strongest electronegative region, dark blue color indicates the strongest electropositive region whereas the green color indicates zero potential.https://doi.org/10.1371/journal.pone.0294511.g012

Fig 13 .
Fig 13.RMSD of the Cα atoms of ROCK2 and the best selected ligand A1 (A) and A2 (B) overtime.(C) RMSD of apoprotein.The left Y-axis shows variation of protein RMSD (blue line represent results) while right Y-axis show ligand variation (red lines) throughout the simulation.https://doi.org/10.1371/journal.pone.0294511.g013

Fig 16 .
Fig 16.Conformational analysis of (A) A1 and (B) A2 systems.Overlays for both protein-ligand complexes at initial (0 ns) represented in blue, middle (100 ns) in yellow and final (200 ns) with grey color frames to analyze conformational changes.https://doi.org/10.1371/journal.pone.0294511.g016 Table in S1 File.Missing residues were added, and model with lowest zDOPE score (-1.16) and highest GA341 score (1.000) was chosen as the best one.